# This code allows to replicate Figures 1-2 in the paper

# It uses an R function for making stacked coefficient plots from different regression estimates
# This function was coded by Piero Stanig, based in part on Yu-Sung Su's coefplot code

stackplot<- function(fit1,
                     longnames=NULL, xlim=c(-.5,0), axis.breaks=NULL,
                     x.label="", y.label="", main.label="", sub.label="",
                     varnames=TRUE,  one.sd=FALSE,
                     lwd.sd=1,max.n=15, offset=.1, smallish=.02,
                     box=FALSE, symmetric=TRUE, cex.var=1.2, cex.pts=1.5,  suppress.intercept=TRUE, 
                     restrict=NULL, suppress=NULL, ind.labels=NULL, pch.type=21, graph_with_custom_axis = NULL
){
  
  # fit: a matrix  
  #		or data frame of coefficients and standard errors
  # longnames: long variable names (vector). If not provided, 
  # 	the variable names from the regression 
  #		that has more predictors are used.
  #restrict: display only some of the coefficients,
  #suppress: do not display dots for a given model  for now, functionality restricted to
  ## matrix input?
  # x.label: label for x-axis
  # y.label: label for y-axis
  # one.sd: if TRUE, plot only the segment for mean plus/minus one s.d.
  # lwd.sd: thickness of the standard deviation bars. 
  # main.label: label for the whole plot
  # sub.label: subtitle. Ignored if estimators=TRUE (see below). 
  # varnames: if FALSE will not plot variable names, default=TRUE
  # box: draws a line between groups of estimates, default=FALSE
  # symmetric: if TRUE makes the plot symmetric around 0, default=TRUE
  # cex.var: the fontsize of the varible names, default=1.2
  # cex.pts: the size of data points, default=1.5  
  
  
  
  
  
  
  if (class(fit1)[1]=="matrix"|class(fit1)[1]=="data.frame"){
    
    coeffs<-as.matrix(fit1)
    
    n.models<-ncol(coeffs)/2
    
    estimators<-FALSE
    
    
    if (is.null(longnames)){
      
      if (!is.null(row.names(fit1))){
        
        longnames<-row.names(fit1)
        
      }else longnames<-" "
      
    }
    
    n.vars<-nrow(coeffs)
    
    max.length<-n.vars
    
    if (max.length>40) {
      
      stop(message = "Can display up to 40 coefficients per model in a plot!\n")
      
    }
    
  }
  
  
  
  
  ####Manage the variable names
  
  if (is.null(longnames)){
    
    longnames<-predictor.names[1:max.length,]
    
    longnames<-t(na.exclude(t(longnames)))[,1]
  }
  
  ####Stack the coefficients and the sd's
  
  coef.mat<-array(NA, c(n.models*n.vars,2))
  
  suppress.mat<-array(NA, c(n.models*n.vars,1))
  
  for (i in 1:n.models){
    
    for (j in 1:n.vars){
      
      coef.mat[((j-1)*n.models)+i,  1] <-coeffs[j,(i*2)-1]
      
      suppress.mat[((j-1)*n.models)+i,  1] <-is.element(i,suppress)
      
      
      coef.mat [i+(j-1)*n.models, 2]<-    coeffs[j,i*2]
      
      
    }
  }
  
  coef<-coef.mat[,1]
  
  
  sd<-coef.mat[,2]
  
  omit.from.plot <- as.numeric(suppress.mat)
  
  
  missing.c<-is.na(coef)
  
  missing.s<-is.na(sd)
  
  coef[missing.c==TRUE]<-0
  
  sd[missing.s==TRUE]<-0
  
  
  
  ###now stack them
  
  n.x <- length(coef)
  
  ####GROUP THE idx SO THAT THE COEFFICIENTS 
  ####ARE CLOSE TO EACH OTHER
  
  idx <- rep(seq(1 , n.x/n.models),each=n.models)
  
  smallsteps<-rep(.4*seq(0,1,length.out=n.models), n.x/n.models) 
  
  idx<-idx+smallsteps
  
  min.x <- round(min(coef-2*sd), 1)-0.1           
  
  max.x <- round(max(coef+2*sd), 1)+0.1                                 
  
  min.y <- round(min(coef-2*sd), 1)-0.3           
  
  max.y <- round(max(coef+2*sd), 1)+0.1                                 
  
  x.scale <- c(min.x+3*offset, max.x+offset)                                            
  
  y.scale <- range(idx)                   
  
  if(offset!=.1){symmetric=FALSE}
  if (symmetric) x.scale<-c(-max(abs(min.x), abs(max.x)), max(abs(min.x), abs(max.x)) )
  
  
  
  subtitle<-""
  
  
  
  if (subtitle=="") subtitle<-sub.label
  
  notshown <- as.numeric(missing.c)+as.numeric(omit.from.plot)
  notshown[notshown==2]<-1
  
  par(mar=c(1+1*as.numeric(estimators==FALSE), 
            2, 2+1*as.numeric(main.label!=""), 
            0))
  
  
  plot(coef, idx, axes=F, type="n",                                     
       
       xlim=xlim, ylim=y.scale+c(-.3, .3), 
       
       xlab=x.label, ylab=y.label,                   
       
       main=main.label, sub=subtitle)                                                        
  
  if(is.null(axis.breaks)){
    if (graph_with_custom_axis == "graph3") {
      
      axis(1, c(-.2, -.1,0))    
    } else {
      axis(1)}
    }
    else{
      axis(1, at=axis.breaks)
    }
  
  abline(v=0, h=0, lty=2, lwd=.5, col="grey40")                                                 
  
  if(n.models==1|length(pch.type)==1){
    
    points(coef, idx-0.25, pch=pch.type, cex=cex.pts*(1-notshown))
    
  } else{
    
    temp.i =1 
    
    for( i.dots  in 1:length(idx)){
      
      points(coef[i.dots], idx[i.dots], pch=pch.type[temp.i], cex=cex.pts*(1-notshown[i.dots]))
      
      temp.i=temp.i+1 
      
      if(temp.i> length(pch.type)){temp.i=1}
      
    }
  }  
  if (box==TRUE){
    
    for (k in 1:n.vars){
      
      segments(x0=-.5, y0=k-.5, x1=0, lty=3)
    }
    
  }
  if (box==TRUE & graph_with_custom_axis == "graph1"){
    
    for (k in 1:n.vars){
      
      segments(x0=.6, y0=k-.5, x1=0, lty=3)
    }
    
  }
  
  if (one.sd==FALSE) segments ((1-notshown)*(coef+2*sd), 
                               idx-0.25, (1-notshown)*(coef-2*sd), idx-0.25, col="Grey80", lwd=lwd.sd)
  
  
  segments ((1-notshown)*(coef+sd), idx-0.25, 
            (1-notshown)*(coef-sd), idx-0.25, lwd=lwd.sd)     
  
  
  
  
  
  
  if (varnames == TRUE) {
      axis(4, .2 + as.numeric(n.models > 1) + (1:n.x), longnames[1:n.x], las = 2, tck = TRUE, pos = -0.515, lty = 0, adj = 1, cex.axis = cex.var)
    
  }
  if (graph_with_custom_axis == "graph1") {axis(4, .2 + as.numeric(n.models > 1) + (1:n.x), longnames[1:n.x], las = 2, tck = TRUE, pos = 0.15, lty = 0, adj = 1, cex.axis = cex.var)
    
  } 
  if (graph_with_custom_axis == "graph3") {axis(4, .2 + as.numeric(n.models > 1) + (1:n.x), longnames[1:n.x], las = 2, tck = TRUE, pos = -0.2, lty = 0, adj = 1, cex.axis = cex.var)
    
  } 
  if(!is.null(ind.labels)){
    
    text(coef, idx-smallish, ind.labels, cex=.8)  
    
  }		
  par(mar=c(5, 4, 4, 2) + 0.1)
}


############################### Figure 1 USA ##############################

cf_elections_us_env = data.frame(outcome= c("Total", "High Income", "Low Income", "China"), 
                                 nocribeta= c(-0.037,-0.077,-0.031, -0.039),
                                 nocrise= c(0.008, 0.009, 0.010, 0.009)
)

cf_elections_us_env <- cf_elections_us_env[nrow(cf_elections_us_env):1,]

varnamez <- c("Total",
              "High Income", 
              "Low Income",
              "China"
)

varnamez <- varnamez[length(varnamez):1]

# Graph will be saved as PDF "Figure_1_USA" in the specified directory

pdf("Figure_1_USA.pdf", width=7, height=9, family="serif")

stackplot(cf_elections_us_env[,c(2,3)], graph_with_custom_axis = "graph3", varnames=FALSE, longnames=varnamez, main = "", xlim=c(-.2,0),
          cex.pts = 0.75, box=TRUE, symmetric = FALSE, cex.var=2.1, offset = 0, pch.type = c(16))

dev.off()

############################### Figure 1 Europe ##############################

cf_elections_eu_env = data.frame(outcome= c("Total", "High Income", "Low Income", "China"), 
                                 nocribeta= c(-0.080,-0.082,-0.212, -0.240),
                                 nocrise= c(0.020, 0.020, 0.088, 0.123)
)

cf_elections_eu_env <- cf_elections_eu_env[nrow(cf_elections_eu_env):1,]

varnamez <- c("Total",
              "High Income", 
              "Low Income",
              "China"
)

varnamez <- varnamez[length(varnamez):1]

# Graph will be saved as PDF "Figure_1_EU" in the specified directory

pdf("Figure_1_EU.pdf", width=7, height=9, family="serif")

stackplot(cf_elections_eu_env[,c(2,3)], graph_with_custom_axis = "graph2", longnames=varnamez, main = "",
          cex.pts = 0.75, box=TRUE, symmetric = FALSE, cex.var=2.1, offset = 0, pch.type = c(16))

dev.off()


############################### Figure 2 USA ##############################

cf_elections_us_att = data.frame(outcome= c("Serious", "Worry", "Environmentalist", "Env_Pri", "More_Jobs"), 
                                 nocribeta= c(-0.028,-0.029,-0.260, -0.005, -0.007),
                                 nocrise= c(0.013, 0.012, 0.117, 0.002, 0.003)
)

cf_elections_us_att <- cf_elections_us_att[nrow(cf_elections_us_att):1,]

varnamez <- c("Serious",
              "Worry", 
              "Environmentalist",
              "Environment Priority",
              "More Than Jobs"
)

varnamez <- varnamez[length(varnamez):1]

# Graph will be saved as PDF "Figure_2_USA" in the specified directory

pdf("Figure_2_USA.pdf", width=7, height=9, family="serif")

stackplot(cf_elections_us_att[,c(2,3)], graph_with_custom_axis = "graph2", longnames=varnamez, main = "",
          cex.pts = 0.55, box=TRUE, symmetric = FALSE, cex.var=2.1, offset = -0.02, pch.type = c(16))

dev.off()

############################### Figure 2 Europe ##############################

cf_elections_eu_att = data.frame(outcome= c("Serious", "Boost economy", "Income insufficient", "Economy gets worse", "Environment priority"), 
                                 nocribeta= c(-0.016,-0.031, 0.264, 0.204, -0.058),
                                 nocrise= c(0.006, 0.006, 0.140, 0.098, 0.032)
)

cf_elections_eu_att <- cf_elections_eu_att[nrow(cf_elections_eu_att):1, ]

varnamez <- c("Serious",
              "Boost Economy", 
              "Income Insufficient",
              "Economy Gets Worse",
              "Environment Priority"
)

varnamez <- varnamez[length(varnamez):1]

# Graph will be saved as PDF "Figure_2_EU" in the specified directory

pdf("Figure_2_EU.pdf", width=7, height=9, family="serif")

stackplot(cf_elections_eu_att[,c(2,3)], graph_with_custom_axis = "graph1", varnames=FALSE, longnames=varnamez, main.label = "", xlim=c(-.1,0.6),
          cex.pts = 0.55, box=TRUE, symmetric = FALSE, cex.var=2.1, offset = 0, pch.type = c(16))

dev.off()


